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ABSTRACT 

Using the convenient second-order interval properties of 
a two-state semi-Markov model for a univariate point proc- 
ess, an automated technique for the estimation of the param- 
eters in the model was researched and discussed. The power 
spectral density of intervals was estimated by the period- 
ogram and a Kolmogorov-Smirnov test of fit was conducted. 

The asymtotic exponential distribution and independence of 
the periodogram points were used to calculate an approximate 
likelihood function. A system of equations was then formed 
to find the maximum likelihood estimates of the parameters. 
Since closed-form solutions for the estimates could not be 
found, an iterative method to stabilize initial guesses of 
the parameter values was attempted with only limited success. 
Results on using Kolmogorov-Smirnov type statistics and the 
spectrum of intervals to test the fit of stochastic process 
models to data have also been obtained. 
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I. INTRODUCTION 



It is in the nature of the Operations Research approach 
to the study of problems to attempt the construction of a 
mathematical model for the problem. Subclasses of mathemat- 
ical models include stochastic, i.e. utilizing random vari- 
ables, and deterministic models. If a stochastic model 
seems appropriate and a general model is proposed, it re- 
mains necessary to estimate parameters of the model from 
data. Parameter estimates, as well as the general form of 
the model, usually come from detailed analysis of data ob- 
served from the problem or process under investigation. 

Several techniques utilizing observed data exist for the 
estimation of parameters for stochastic models. Typically 
the methods of moments or maximum likelihood are used and 
usually yield estimates with some desirable properties. 
Methods such as these frequently require the simultaneous 
solutions to a system of equations in order to find esti- 
mates. A number of computer approximation routines have 
been developed for the solution of such systems, but their 
usefulness seems limited. 

One proposed stochastic model provided the impetus for 
this research. Lewis and Shedler [1973], while studying 
page reference patterns in a demand paged computer system, 
formulated a univariate two-state semi-Markov model for the 
process of page exceptions. Page exceptions occur because a 
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computer program which is in execution has been stored in 
blocks of storage called pages. Some of these pages must be 
in core storage for the program to be executing, while the 
remaining pages may be located on peripheral storage de- 
vices. Following the execution of each instruction a page 
is referenced which contains the next instruction. If this 
referenced page is in core storage execution continues; how- 
ever, if the referenced page is not in core storage then 
execution is interrupted and the referenced page must be 
read into core storage. This type of interruption is refer- 
red to as a page exception. Data for this process was gen- 
erated by counting the number of page references occurring 
between page exceptions. Lewis and Shedler [1973] discussed 
their procedure for estimating parameters which they de- 
scribed as an ad hoc method, and concluded that there was a 
need to formalize the parameter estimation procedure. 

The purpose of the research in this thesis was to uti- 
lize the convenient second-order interval properties of a 
univariate two-state semi-Markov process to produce an auto- 
mated, computer programed, technique for the estimation of 
parameters for the model. This was desirable because the ad 
hoc method used by Lewis and Shedler [1973] was very time- 
consuming and there exists a considerable body of page ex- 
ception data which it is desired to analyze. The basic 
procedure was to calculate an estimate for the power spec- 
tral density of the process, namely the periodogram, and 
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utilize an approximate method of maximum likelihood to esti- 
mate the parameters. 

It will be seen that the proposed procedure did not work 
as well as hoped, but the problems which arose pointed up 
other possible attacks on the problem. It should also be 
noted that model fitting and parameter estimation for these 
point processes is almost a completely open field. 
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II. BACKGROUND ANALYSIS 



A. TWO-STATE SEMI-MARKOV MODEL FOR UNIVARIATE POINT PROCESS 

Excellent discussions of this model can be found in Cox 
and Lewis [1966, Ch. 7] and Lewis and Shedler [1973]. Those 
discussions are summarized here for continuity of exposition. 

Let the sequence of random variables { X^ , i=l, . . . ,N} be 
interevent times, i.e. is the interevent time between 
event (i-1) and event (i) . In order that a discussion of 
equilibrium distributions may be avoided it was assumed that 
a hypothetical event has occurred at time zero, so that X : , 
the interval between time zero and the first event, is an 
observation from the same process as the remainder of the 
sequence, i.e. there is no length-biased sampling [Cox and 
Lewis 1966,Ch.4] included. 

Now suppose there are two types of intervals but that 
the interval type is not observable, i.e. a univariate point 
process. The two interval types have probability mass func- 
tions (p.m.f.) Pi (x) and p 2 (x) , respectively, with transi- 
tions between types described by a two-state Markov chain 
with matrix 




That is, given that X^ has p.m.f. p 2 (x) then X^ +i has p.m.f. 
p 2 (x) with probability a 2 and p.m.f. p } (x) with probability 
l-a 2 , independent of the history of previous intervals, etc. 
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The vector of steady-state probabilities tt = ( tt x tt 2 ) 
associated with the transition matrix A results from the 
solution of the matrix equation tt = ttA and it follows that 



1-a 



1-a 



TTi = 



tt 2 = 



1 



2-ai-a2 2-aj-a 2 

If and y 2 ,0 2 2 are t ^ le mean and variance for intervals 

with p.m.f. pj(x)'and p 2 (x) , respectively, the steady-state 
marginal results for intervals between events in the univar- 
iate process, i.e. interval type not known, are as follows: 

P (x) = TTjPjfx) + ir 2 p 2 (x) , 

y = E (X) = iTiyii + ir 2 y 2 , 

a 2 = var(X) = tt j 0 j 2 + ir 2 a 2 2 + tt 1 tt 2 ( y 1 — y 2 ) 2 . 

The serial correlation coefficients of lag k, p^, for the 
intervals are of the form m B^/ 0 2 where, for k=l,2,..., 
m = tt 1 it 2 ( y ! — y 2 ) 2 B = a x + a 2 -1 . 

From these coefficients the positive portion of the power 
spectral density may be computed, 



P (co ) = tt (1 + 2 E p.cos kco_) . 



+ n 



k=l 



n 



The closed-form solution to the infinite series is given by 

Jolley [1961, series #545] yielding 

0 2 mB (cos w ) - B 

P+(%) = — [1 + 2 — { }]. (1) 

TT 0 2 1+B 2 -2BcOS 0) n 

The beneficial feature of the power spectral density for 
this model is that it only depends upon the mean and vari- 
ance of each of the two probability distributions and not 
on the complete distributions, and is thus fairly robust. 
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The count spectrum [Cox and Lewis 1966, Ch.4] on the other 
hand, depends on the complete distributions. 

B. PARAMETER ESTIMATION TECHNIQUE 

Lewis and Shedler [1973] used a modified method of mo- 
ments approach in order to estimate the parameters in their 
model. The standard method of moments procedure for param- 
eter estimation is to calculate theoretical moments in terms 
of the unknown parameters and equate them to empirical ob- 
servations of these moments. An alternative to this method 
is the method of maximum likelihood. In this method param- 
eter values are selected which maximize the joint probabil- 
ity density of the observed data. To accomplish this there 
is a need for some distributional assumptions. However, 
even for a simple model such as the univariate two-state 
semi-Markov model discussed here it is not possible to write 
down the joint density of the intervals. Thus the following 
approximate technique was proposed and tried. 

It is known, Cox and Lewis [1966, Ch. 5], that an estimate 

of the power spectral density P (oa n ) at w , the periodogram 

I (n) , is in general asymtotically exponentially distributed 

[Olshen 1967] . The periodogram is an unbiased estimate, 

i.e. E[I(n)]=P(co ); however, it is not a consistent estimate 

n 

since the variance of the exponential distribution is equal 
to the square of the mean, i.e. the variance does not de- 
crease with increased sample sizes. Moreover, for n t not 
equal to n 2 the periodogram points I(nj) and I(n 2 ) are 
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asymtotically independent. Thus for finite sample size N an 
approximate likelihood function may be written by assuming 
the periodogram points are independent with exponential dis- 
tributions having mean value P(oo n ). This is the technique 
explored in this thesis. 

The definition and development of the periodogram re- 
quires the finite- Fourier transform. 

The finite Fourier transform was discussed by Cooley, 
Lewis and Welch [to be published in 1974]. Let {Y(j), 
j=0, . . . ,N-1} be a sequence of N real numbers. The finite 

Fourier transform of Y(j) is then 
N-l _ 2irin j 

a(n) = i E Y(j)e N , n=0,...,N-l. 

N j=0 

This sequence of complex numbers may also be written in the 
form 

N-l N-l 

a(n) = — E Y ( j ) cos ( joi ) -i E Y(j)sin(jio ) , 

N j=0 n j=0 

where to = 27m, n=0 , 1, . . . ,N-1 . The periodogram I (n) is then 
n 

N | a (n) | 2 

I (n) = 2 tt , n=0, . . . , 



or 

N-l N-l 

I (n) = { E Y(j)cos(jo) ) } 2 + { E Y(j)sin(jto )} 2 

j=0 _ n 

2ttN 

The periodogram is an even, periodic, function and hence has 
only [N/2J+1 distinct values, where [N/2] is the integer 
part of N/2. Hereafter in the discussion N will refer to an 
even integer. Let I + (n)=2l(n) be the estimate for P + (co n ). 
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It is easily seen that I + (0) is proportional to the 
square of the arithmetic average of the observed data; thus 
no new information is obtained from I + (0). Since N/2 is an 
integer w N /2 =7r an< ^ I + (N/2) is proportional to the square of 
an alternating summation of the data. Both I + (0) and 
I + (N/2) were ignored in what follows, thus leaving (N/2)-l 
periodogram points. It should be added that these two peri- 
odogram points, suitably normalized, have asymtotic x 2 dis- 
tributions with one degree of freedom and not an exponential 
distribution. 

Now there is sufficient information to begin the approx- 
imate maximum likelihood search for parameter estimates. 

The parameters of this model that need estimation are the 
mean and variance of each marginal distribution and the two 
transition probabilities aj and a 2 - As a vector these pa- 
rameters will be labeled £=(y 1 ,a 1 2 ,y 2 ' a 2 2 ' a i' a 2 ) an d indi- 
vidually, to simplify notation, as 0 j , j=l , 2 , 3 , 4 , 5 , 6 , to 
stand for the parameter as an element of the vector 0_. 

The approximate likelihood function can be written as 

(N/2 ) -1 1 -I + (n)/P.(m ; 6) 

L ( 0 ) = n PTTw ; 0) e 

n=l n “ 



which is equivalent to 



L(0) 



/ (N/2 ) -1 

n p 
V n=l 






fN/ 2 )" 1 1+ (n) 

e n=1 • 



A more simple function to work with, which has the same max- 
imum as L(0_), is’ the log likelihood function LL(0^=ln L (0_) , 
where In symbolizes the natural logarithm. 
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Then 

(N/2 ) -1 (N/2) -1 

LL ( 0 ) = - E ln{P.(w ; 6 ) ) - E I +^ n > 

n=l + n n=l P + (u n ;6) . 

In the typical mathematical approach to finding an un- 
constrained maximum of a function, it is a necessary condi- 
tion that all of the first-order partial derivatives of the 
function, with respect to the unknown parameters, be equal 
to zero, i.e. that 

0 = 3LL (1) = LL- = E * J + (h) ~ P + P . , j=l , . . . , 6 (2) 

36 j n=l P + (w n ; 0_) J 

where Pj and LLj are the first-order partial derivatives of 
P + (w n ?£) and LL (B) , respectively, with respect to parameter 
0j. This process results in six equations and six unknown 
parameters. Parameter estimates are found by simultaneously 
solving the system of equations for each of the parameters, 
although this may not yield a unique maximum. If the system 
is of a simple form it may be possible to get at least a few 
closed-form solutions which will reduce the size of the sys- 
tem. 

Once the parameter estimates have been found it is nec- 
essary to show that a maximum has been achieved. A suffi- 
cient condition for a maximum is that the matrix of second- 
order partial derivatives be negative definite. The final 
phase in this approximate likelihood estimation process is 
to verify the predictability of the model. The verification 
may be done, using the estimated parameters, by calculating 
other theoretical properties of the model, such as the spec- 
trum of counts discussed by Cox and Lewis [1966, Ch.4], 
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which may then be compared with the corresponding empirical 
properties of the data. Note that the utility of the spec- 
trum of intervals in the approximate likelihood analysis is 
that it does not depend on the complete distributional form 
for pi (x) and p 2 (x) while the spectrum of counts does. It 
will be seen later, however, that this independence leads to 
ill-conditioning in the solution of the maximum likelihood 
equations . 
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III. EXPERIMENTAL APPROACH 



The original data, analyzed by Lewis and Shedler [1973] , 
was not available for this research. In view of this fact 
and since the purpose of the research was to evaluate the 
effectiveness of the previously described technique for pa- 
rameter estimation it was felt that data observed from a 
model with known parameters would better aid the evaluation 
process. With this in mind a simulation of the model de- 
scribed by Lewis and Shedler [1973] was constructed for the 
purpose of generating such data. 

A. UNIVARIATE TWO-STATE SEMI-MARKOV SIMULATION 

The simulation, as well as the model, was subdivided 
into three major subsections. The state transition matrix A 
was one subsection and the two distributions for intervals 
were the remaining two subsections. Lewis and Shedler [1973] 
postulated a geometric distribution for the long intervals 
and a negative binomial distribution for the shorter inter- 
vals. The parameters used for the simulation were those 
calculated by Lewis and Shedler [1973] . 

A Monte Carlo simulation, such as this, required a 
pseudo-random number generator with favorable serial corre- 
lation properties. Learmonth and Lewis [1973] discussed 
such a generator called SRAND. SRAND returns an observation 
from a standardized uniform distribution on the interval 
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(0,1) . SRAND is a multiplicative generator with a multi- 
plier of (7 s ) and a modulus of (2 31 -1). 

The geometric distribution is of the form 

Pj(x) = p^" 1 (1-pj), 0<p 1 <l; x= 1,2,..., 
with a mean y^l/d-pj) and variance c x 2 =p 1 / (1-Pj ) 2 . Uti- 
lizing the survivor function of the geometric distribution, 
i.e. prob{X>x}=p 1 x , x=l,2,... , a generator of geometric 
variates was obtained. It was of the form 
x = T {In (R) /In (pj) } , 

where R was an observation from SRAND and the symbol f {b} 
signified the smallest integer greater than or equal to b. 

The negative binomial distribution is of the form 



P 2 (*) 




d-P 2 ) k , 



0<p 2 <l, k>0, x=l,2,... , with mean y 2 =l+{kp 2 / (l-p 2 ) } and 

variance o 2 2 =kp 2 / (l-p 2 ) 2 . Let x|x, denoting X given a fixed 
value of X, be distributed as a Poisson random variable with 
parameter X. Now let X have a gamma distribution with pa- 
rameters k and ri/ 



f (X) 



n k x k ~‘ e-^ 
r(k) 



k>0; X>0; 



ri > 0 . 



It can be shown using generating functions that the uncondi- 
tioned X has a negative binomial distribution with parame- 
ters k and p 2 =l/(l+n)* 

To calculate a gamma variate with a parameter k, a posi- 
tive real number, it was necessary to employ Johnk's tech- 
nique [1964] for generating variates with the fractional 
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part of k. Let k be the integer part of k, if k>l, or zero 
if k<l, and let k be the fractional part of k. The sum, X x , 
of k exponentially distributed random variables with param- 
eter n has a gamma distribution with parameters k and n* In 
Johnk's technique let U lf U 2 and U 3 be independent and iden- 
tically distributed observations from a uniform distribution 
on interval (0,1). such that 

Y = U/A + u 2 ly/(1- ^ <1 . 

If Y>1 new observations for U x and U 2 should be obtained. 
Then for Z=U i 1 ^/Y and E=-ln U 3 , X 2 = ( Z X E ) /n has a gamma dis- 
tribution with parameters k and g. Finally X=X x +X 2 has the 
required gamma distribution with parameters k and n* 

The generation of Poisson random variates with parameter 
X was accomplished by letting X be equal to the largest in- 
teger n such that, for a sequence of independent identically 
distributed uniform random variates (U^) from the interval 
( 0 , 1 ), 

U.xU X...XU >e -X . 
n 

If Uj<e - ^ then X=0. X is then distributed as a Poisson var- 
iate with parameter X. 

B. _ CALCULATION OF THE PERIODOGRAM 

The finite Fourier transform discussed in section II. B 
above requires on the order of N 2 complex operation pairs, 
i.e. a multiplication and an addition. For large N this 
can be very costly in terms of calculation time. Cooley, 
Lewis and Welch [1970] discussed the use of a fast Fourier 
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transform algorithm which only requires on the order of 

N(r!+...+r ro ) complex operation pairs where N= (r t x. . . x r m ) , 

i.e. r is a factor of N. The International Mathematical 
m 

and Statistical Library [1973 revision] contains a computer 
subroutine, FFTR, which computes the fast Fourier transform 
of a real data sequence. For N=820, as in this research, 
the fast Fourier transform algorithm used only six percent 
of the number of complex operation pairs required by the 
straight-forward calculation method. Thus a significant 
savings in computer operating time was realized. 

Utilizing previously described equations the periodogram 
I + (n) was computed and then used in a test of fit to the 
power spectral density P + (w n ). Cox and Lewis [1966, Ch.6] 
described a test based on the uniform distribution. While 
the periodogram has, asymtotically , an exponential distri- 
bution with mean P + (w n ), the quantity I + (n) /P + (u n ) has an 
exponential distribution with mean one. This is true for 
each of the (N/2)-l periodogram points. If all (N/2)-l of 
these quantities are summed the total gives an interval of 
length over which there are (N/2)-2 points dispersed. The 
intervals between these points are each, hypothetically, an 
observation from a unit exponential distribution, i.e. the 
points form a Poisson process. It is a well known fact of 
the Poisson process that given M points are in an interval 
the M points are dispersed uniformly over the interval. 
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Thus, the sequence {U^, i=l, . . . , (N/2) -2 } , where 

nSl I+(n)/P + (u n ) 

U (i) — (N/2 ) -1 

„?1 V n)/P + <“„> 

are uniform order statistics. The empirical cumulative dis- 
tribution function for these quantities was then compared 
with the uniform cumulative distribution function using the 
Kolmogorov-Smirnov test. The null hypothesis is that the 
sequence is formed of uniform order statistics, while 

the alternative hypothesis remains unspecified. Lilliefors 
[1969] found that the critical values of the K-S test are 
too conservative when testing using exponential distribu- 
tions where the mean has been estimated, as in this case. 

Too conservative means that the listed critical value for a 
level of significance a has actually a level of significance 
less than a. If the above test, with modified percentage 
points, accepts the null hypothesis then the assumption of a 
semi -Markov model for the data has been justified. 

In order to test the periodogram it was necessary to 
know P + (a) n ) . As discussed earlier the correlation coeffi- 
cient of lag k,Pj., is /a for this model. Let y(0), 

Y(l) and y(2) be estimates of the variance and covariances 
of lags one and two, respectively, for the intervals. Then 
Y (1) = c 2 p a = m3 

and 

y (2) = o 2 p 2 = m3 2 
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Solving simultaneously for in and B, the estimates of m and 3 

are 

3 = y (2 ) and m = y 2 (1) 

ytlT yT2T . 

From (1) an estimate for P + (w n ) was 

1 „ (cos to ) -r 

p + (w n ) = it {y (0) + 2m 3 } 

1+ 3 2_2 3 cos w n 

These estimates were then used in the computations for the 
sequence {U^}. 

C. SOLVING SIMULTANEOUS EQUATIONS FOR THE PARAMETERS 

The system of equations defined by (2) and (1) was ex- 
tremely complex, with no hope of finding a closed- form solu- 
tion for any of the parameters. The system was reduced, 
however, by noting from the geometric distribution assump- 
tion that the variance for the long intervals was a function 
only of the parameter pi, which also was the only parameter 
in the mean. It was a simple matter then to find the vari- 
ance as a function of the mean which then reduced the system 
to only five unknown parameters. The system was still com- 
plex and required some iterative method for solution. 

Rao [1965] suggested an iterative method which he called 
the Method of Scoring. He called LL j the jth efficient 
score. The approach for this method was to assume some ini- 
tial trial solution. Using a first-degree Taylor's expan- 
sion of the efficient scores about the trial solution, a 
system of linear equations was derived from which an addi- 
tive correction to the trial solution was found. The 
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iterations were repeated until the additive corrections be- 
came negligible. 

Specifically, let 0 1 °,...,0 5 ° be the trial values for 
the unknown parameters. From (2) 

3LL ( 0.) 5 3 2 LL (£.) 

LL. * 30-; 0 + .1 (0i-0i°) 30^°30. 0 , j=l,...,5 , 

J J 1=1 J 1 

where 3LL (£) /30 j °=S j 0 the first-order partial derivative of 
LL(0_) with respect to 0j, evaluated at 0 j 0 . Let 60j=(0j-0j°) 
and also let 3 2 LL (0j / ( 30 j 0 30^ ) =Tji 0 . Then 
5 

- Z T..°50 i = S,°, j=l , . . . , 5 
i=l - ,1 J 

was a system of linear equations with five unknowns. In 
matrix notation the system had the form -T<S 9= S . Finally, 
the additive corrections were obtained from the equation 
6 9=-T~ 1 S where T -1 was the inverse of the matrix T, assuming 
T was nonsingular. The new trial solution then became 
0/ = 9 °+6 9 . 

Rao [1965] explained that the variance of the final es- 
timate 9jf of 0j was approximated by the jth diagonal ele- 
ment of the matrix (-T -1 ). Recalling that the matrix of 
second-order partial derivatives of LL (0J should be negative 
definite, then -T and (-T -1 ) were both positive definite. 

In order to apply the method of scoring it was necessary 
to determine initial estimates of the parameters. The mean 
and variance of the intervals and the parameters m and 0 all 
have been estimated. Utilizing the marginal properties of 
the model and the method of moments a system of four 
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equations was developed which was of the form 

X = (l-5 2 )fii + (l-5!)y 2 

0^1 ; 

Y(0) = (l-5 2 ) (Pi 2 -yi) + (l-ai)g 2 2 + m ; 

( 1 - 6 ) 

6 - 5i+a 2 -l ? 

m = (l-5i) (l-5 2 ) (y t -y 2 ) 2 
( 1 - 6 ) 2 

where X was the estimate for the mean of the intervals and 

the quantity (yj 2 -^) was the estimate for Oj 2 after making 

the geometric distribution assumption. From this system 

initial estimates for four parameters, as functions of the 

fifth parameter, were found and had the form 

Xyj - X 2 - m 
Vz = I 

(yj- x ) 

(1-6)X - (y 2 -BiJ 1 ) 

5 i = : — : ; 

(y 1 -U 2 ) 

a 2 = 6+1-cti ; 

(1-6) (Y(O)-m) - (l-a 2 ) (y j 2 -y x ) 



It only remained then to estimate parameter y x . 

Lewis and Shedler [1973] explained that their estimate 
of P! involved an eyeball judgement of where linearity began 
in the tail of the log survivor function of the data. This 
linearity in the tail led to the postulate of the geometric 
distribution for the long intervals. Since only an initial 
estimate was needed, their method of estimating yi was uti- 
lized again. The value of the interval where linearity 



22 



began was subtracted from each greater interval and the 
arithmetic average of these intervals was then taken as the 
estimate of yi. Now all parameters had been initially esti- 
mated and the iterative method of scoring was applied to 
stabilize these estimates. 
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IV. RESULTS AND CONCLUSIONS 



As the research for this thesis progressed two areas 
developed results which need to be discussed. The first of 
these was the test for justification of the exponential dis- 
tribution assumption for the periodogram points. A subrou- 
tine called KSTEST was written to conduct this test. Part 
of the output of this subroutine was the Kolmogorov-Smirnov 
statistic which then was compared to a critical value from 
the distribution proposed by Lilliefors [1969] for the case 
where the mean of the exponential distribution had to be es- 
timated. The 0.99 quantile of that distribution, i.e. a one 
percent level of significance, was 1.25. It was noted that, 
at this level, of four thousand trials made approximately 
six percent were rejected as not having produced periodo- 
grams from a semi-Markov model. 

In addition to testing the hypothesis for each simula- 
tion another benefit was received. Since the testing did 
not strictly conform to that discussed by Lilliefors [1969] , 
because each periodogram point had a different mean, it was 
felt that, for this case, quantiles of the distribution 
should be estimated. The four thousand data points of the 
statistic were obtained from four computer runs each contain- 
ing one thousand simulations. For each run, the data was 
divided into ten sections in serial order, i.e. the first 
one hundred points were the first section, etc. The ele- 
ments of each section were ordered and the 0.80, 0.85, 0.90, 
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0.95 and 0.99 empirical quantiles were observed. This re- 
sulted in ten observations for each quantile from which a 
mean and variance were estimated. Lastly, the entire data 
for the run was ordered and the five quantiles were observed. 
Thus, for each run, each of the five empirical quantiles had 
been observed and had an estimated mean and variance. Fin- 
ally a mean of the four overall observations for each quan- 
tile and a mean of the section means were computed. The 
results are shown in Table I. 

Lilliefors [1967] discussed the Kolmogorov-Smirnov test 
for normal data and calculated numerically the quantiles of 
the test statistic for the case where the mean and variance 
of the normal distribution must be estimated. These quan- 
tiles are included for comparison. 

The second of the two significant areas was the estima- 
tion of parameters. The subroutine ESTIM8 was written, in 
double precision, to utilize the method of scoring for pa- 
rameter estimation. As a result of the use of the subrou- 
tine several potential hazards to the proposed technique 
became visible. 

The first of these hazards was the disparity between 
magnitudes of the five unknown parameters. Three of these 
are means and variances and the other two are probabilities, 
which are always less than or equal to one. This problem 
became apparent when the magnitudes of the scores and the 
elements in the matrix of second-order partial derivatives 
were seen. An attempt to correct this problem was made by 
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TABLE I 



Source Level of Significance 





0.20 


0.15 


0.10 


0.05 


0.01 


Usual quantiles 


1.07 


1.14 


1.22 


1.36 


1.63 


Lilliefors quantiles 


0.86 


0.91 


0.96 


1.06 


1.25 


Run 1 


0.74 


0.79 


0.96 


1.66 


8.56 


Mean 


0.73 


0.80 


0.92 


1.59 


6.04 


Variance 


0.002 


0.005 


0.01 


0.31 


13.24 


Run 2 


0.75 


0.85 


1.03 


2.05 


6.66 


Mean 


0.75 


0.85 


1.03 


1.86 


7.11 


Variance 


0.003 


0.01 


0.04 


0.45 


25.35 


Run 3 


0.74 


0.80 


0.90 


1.31 


6.79 


Mean 


0.73 


0.79 


0.91 


1.80 


4.90 


Variance 


0.001 


0.002 


0.01 


1.31 


9.87 


Run 4 


0.73 


0.82 


0.96 


1.72 


8.45 


Mean 


0.74 


0.83 


0.95 


1.64 


5.90 


Variance 


0.01 


0.02 


0.03 


0.62 


12.51 


Mean of Runs 


0.74 


0.82 


0.96 


1.69 


7.62 


Variance of Runs 


0.00 


0.00 


0.001 


0.03 


0.27 


Mean of Means 


0.74 


0.82 


0.95 


1.72 


5.99 


Variance of Means 


0.00 


0.00 


0.001 


0.01 


0.21 


Lilliefors normal 












quantiles 


0.736 


0.768 


0.805 


0.886 


1.031 
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dividing the three large parameters by the overall mean of 
the intervals. This also favorably affected the partial 
derivatives involving these parameters. The desired effect 
was achieved in that the gap between magnitudes was nar- 
rowed; however, the parameter estimates that resulted from 
this modification were only about one percent different from 
the parameters achieved earlier, so apparently the disparity 
created no significant problem. 

The second of these problems was that the final param- 
eter estimates, overall, appeared to have little relation- 
ship to the marginal parameter values from which the data 
was generated. Similarly, parameter estimates for two sets 
of data differed greatly in magnitude and at times in sign, 
even when the periodogram was accepted as a close fit to the 
power spectral density. Differences in sign were extremely 
disturbing since all of the parameters were expected to be 
greater than zero. 

A third problem, related to the second, was that the re- 
sults failed, numerically, to establish that the matrix of 
second-order partial derivatives was negative definite. 
Similarly the negative inverse of that matrix could not be 
shown to be positive definite. This problem indicated that 
either a maximum had not been achieved, even though a cut- 
off criterion of 10 -10 was used to test for convergence, or 
that due to round-off error the properties of a maximum 
could not be detected. With a smaller cut-off criterion the 
process would not converge and had to be terminated. 
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Finally, in a few instances when four of the five un- 
known parameters appeared close to the simulation parameters 
the initial value of the fifth parameter was changed and the 
subroutine was restarted. The parameters would again con- 
verge; however, the final values in some cases changed dras- 
tically, even to the point of changing sign. 

Some of these • problems may have been caused by an ill- 
conditioned system of equations, while others might be due 
to the lack of a powerful iterative technique for the solu- 
tion of a system of equations that has, perhaps, poor initial 
estimates. In any case it should be clear that the use of 
second-order properties of a model might simplify or at 
least aid the parameter estimation process. One proposed 
modification to the technique discussed in this thesis was 
to use a mixture of the method of moments approach on the 
marginal distribution of the intervals and the maximum like- 
lihood approach on the second-order properties to estimate 
parameters. 

In conclusion it should be recalled that model fitting 
and parameter estimation for univariate point processes is 
almost a completely open field and that attempts, even un- 
successful ones, are needed in order to break-through the 
barrier of inadequate methodology. 
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COMPUTER SUBPROGRAMS 



SUBROUTINE MODEL 

PURPOSE: 

THIS SUBROUTINE CONTROLS THE SIMULATION OF A 
UNIVARIATE TWO-STATE SEM I -MARKCV POINT PRGCESS. 

THE STATE ONE INTERVALS HAVE A GEOMETRIC 
DISTRIBUTION AND THE STATE TWO INTERVALS HAVE 
A NEGATIVE BINOMIAL DISTRIBUTION. 

USAGE: 

1. ARGUMENTS: 

X - OUTPUT VECTOR OF INTEREVENT TIMES (REAL*8) 

SIZE - INPUT LENGTH OF VECTOR X (INTEGER) 

IX - INPUT RANDOM NUMBER SEED (INTEGER) 

DIST1M - INPUT MEAN OF THE GEOMETRIC 
DISTRIBUTION (REAL*8) 

DIST2M - INPUT MEAN OF THE NEGATIVE BINOMIAL 
DISTRIBUTION (REAL*8) 

DIST2K - INPUT PARAMETER K IN THE NEGATIVE 
BINOMIAL DISTRIBUTION (REAL*8) 

A 1 - INPUT PARAMETER IN THE TRANSITION 
MATRIX (REAL*8) 

A2 - INPUT PARAMETER IN THE TRANSITION 
MATRIX (REAL*8) 

2. REQUIRED SUBPROGRAMS: 

INTEGER FUNCTION GEOMET 
INTEGEP FUNCTION NEGBIN 
SUBROUTINE OVFLOW (NPS ROUTINE) 

SUBROUTINE SR AND (NPS ROUTINE) 

SUBROUTINE SNORM (NPS ROUTINE) 

SUBROUTINE SEXPON (NPS ROUTINE) 



SUBROUTINE MODEL ( X , S I ZE , I X , D I ST 1M , D I ST2M , D I ST2K, A 1 , A2 ) 
IMPLICIT REAL *8 (A-H,K,0-Z) 

REAL** R 

INTEGERS GEOMET, STATE.SIZE 
DIMENSION ALPHA( 2 ) ,X( S IZE) 

COMMON ETA , K, PARAMG 
CALL OVFLOW 

PARAMETERS FOR NEGBIN 
K = D I ST 2 K 

ETA=K/( DIST2M-1 .ODO) 

PARAMETER FOR GEOMET 

PARAMG=DLOG( 1 .ODO-1 .ODO/DIST1M ) 

STATE ONE TRANSITION PROBABILITIES FOR MATRIX 
ALPHA ( 1 ) = A I 
ALPHA ( 2 )= I .ODO-A2 

SELECT INITIAL STATE FROM STEADY-STATE PROBABILITIES 
S TAT E= 1 

CALL S R AND ( I X , R , 1 ) 
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IF(R.LE.((1.0D0-ALPHA(1))/(1.0D0-ALPHA(1)+ALPHA(2) ) 
O) ST AT E=2 

COMPUTE 'SIZE' INTEREVENT TIMES 
DO 2 1*1, SIZE 

ENTER MATRIX AND DETERMINE TYPE OF NEXT INTERVAL 
CALL SRANDI IX, R , 1 ) 

IF(DBLE(R) .LE.ALPHA(STATE)) GO TO 1 

PICK TYPE TWO VARIATE 
STATE=2 

X(I )=DFLOAT(NEGBIN(IX) ) 

GG TO 2 

PICK TYPE ONE VARIATE 

1 STATE=1 

X ( I )=DFLOAT(GEOMET (IX) ) 

2 CONTINUE 
RETURN 
END 



INTEGER FUNCTION GEOMET 

A. .PURPOSE: 

THIS FUNCTION GENERATES VARIATES FROM THE 
GEOMETRIC DISTRIBUTION WHICH IS CF THE FORM 

X-l 

M(X)={ 1-PJ*( P) ;0<P<l;X=l,2»... 

B. USAGE: 

THIS FUNCTION WAS WRITTEN TO BE USED WITH 
SUBROUTINE MODEL. 

P=l-l/( 1-DIST1M) 

PARAMG=DL0G(P) 



INTEGER FUNCTION GEOMET ^4 (IX) 

IMPLICIT R EAL*8 (A-H,K,0-Z) 

REAL*4 R 

COMMON ET A t K t PARAMG 

CALCULATE VARIATE 
CALL S RAND ( I X * R » 1 ) 

RATIO = DLOG(DBLE( R > I /PARAMG 
IRATIO = IDINT(RATIO) 

ROUND UP IF NON-INTEGER 

IF (RATIO-DFLOAT( IRATIO). GT.O.ODO) IR AT I 0= I RAT 10+ 1 

8 RETURN VARIATE 
G EGMET = IRAT I 0 
RETURN 
END 
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INTEGER FUNCTION NEGBIN 



A. 



6 . 



PURPOSE: 

THIS FUNCTION GENERATES VARIATES FROM THE NEGATIVE 
BINOMIAL DISTRIBUTION WHICH IS OF THE FORM 



(K-X-2) 
M ( X ) = ( ) 

( X-l ) 



K X-l 
*(1-P) *(P) 



;K>0;0<P<l;X=l,2, . .. 



USAGE: 

THIS FUNCTION WAS WRITTEN TO BE USED WITH 
SUBROUTINE MODEL. 

P=l/ ( 1 + ETA ) 



INTEGER FUNCTION NEGB IN*4. ( IX) 

IMPLICIT R EAL *8 (A-H,K,G-Z) 

REAL*4 U( 1 00 ) , E » NORM 
COMMON ETA » K» PARAMG 
G AMMAD= 0 . 0D0 
GAMMA I = 0 . 0 DO 

DETERMINE INTEGER AND DECIMAL PARTS OF K 
I K=IDI NT ( K ) 

DK=K-DF L0AT( IK) 

CALCULATE, IF REQUIRED, GAMMA IK VARIATE 
FROM SUM CF IK UNIT EXPONENTIAL VARIATES 
IF(K.LT.l.ODO) GO TO 9 
E T=0 . 0D0 
DO 8 M=l, IK 
CALL SEXPON (IX,E,1) 

ET=ET+DBLE(E) 

8 CONTINUE 
GAMMAI =ET/ETA 

CALCULATE, IF REQUIRED, GAMMA DK VARIATE 
USING JOHNK ' S METHOD 

9 IF(DK.LE.O.ODO) GO TO 11 

10 CALL S RAND (IX, U, 2) 

UK1 = D3LE(U(1) )**( 1.0D0/DK) 

UK2=DBLE(U (2 ) )#*(1.0D0/(1 .ODO-DK ) ) 

ZZ=UK1+UK2 

IF(ZZ.GE.l.ODO) GO TO 10 
CALL SEXPON (IX, E, 1) 

GAMMA D= (UK1/ZZ)*DBLE(E)/ETA 

TOTAL GAMMA VARIATE 

11 GAMMA=GAMMAD+GAMMAI 
I GAMMA= ID I NT (GAMMA) 

IF( IGAMMA.GE.100) GO TO 50 

CALCULATE POISSON VARIATE, DIRECTLY 
NN = 0 

UT = 1 . 0 DO 

EGAMMA=DEXP (-GAMMA) 

20 CALL S RAND ( IX, U, 100) 

DO 30 M=1 , 100 
UT=UT#DBLE(U(M) ) 

IF(UT.LE.EGAMMA) GO TO 40 
30 CONTINUE 
NN=NN+ 1 00 
GO TO 20 
40 N = M- 1 + N N 
GO TO 60 

CALCULATE POISSON VARIATE, USING NORMAL APPROXIMATION 
50 CALL SNORM (IX, NORM, 1) 
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N = IDI NT (GAMMA-0.25D0+( NORM/ 2 . ODO ) **2+N0 RM* 
C ( OSQRT ( GAMMA+O . 125D0 ) ) ) 

RETURN VARIATE 
60 NEGBI N= N+l 
RETURN 
END 



SUBROUTINE KSTEST 

A. PURPOSE: 

THIS SUBROUTINE CALCULATES THE PERIODOGRAM OF 
INTERVAL DATA, ESTIMATES THE POWER SPECTRAL 
DENSITY PSD AND TESTS THE FIT OF THE PERI CDOGRAM 
TO THE PSD. 

E. USAGE: 

1. ARGUMENTS: 

X - INPUT VECTOR GF INTERVAL DATA (REALMS) 

I VEC - OUTPUT VECTOR OF DESIRED PERIODOGRAM 
POINTS (REAL*8) 

MEAN - OUTPUT MEAN OF INTERVAL DATA (REAL*8) 

VAR I AN - OUTPUT VARIANCE OF INTERVALS (REAL*3) 

XKSDN - OUTPUT KOLMOGOROV-SMI RNOV STATISTIC 
FOR TEST OF FIT (REALMS) 

SIZE - INPUT LENGTH OF VECTOR X (INTEGER) 

MHAT - OUTPUT ESTIMATE OF COVARIANCE PARAMETER 
M (REAL*8) 

BHAT - PUTPUT ESTIMATE OF COVARIANCE PARAMETER 
BETA (REALMS) 

2. REQUIRED SUBROUTINE: FFTR (IMSL ROUTINE) 

3. CAUTION: VECTOR X OF INTERVALS IS DESTROYED 

BY FFTR. 



SUBROUTINE KSTEST ( X , I V EC , ME AN , V AR I AN, KS DN , S IZ E , MHA T , 
C BHAT ) 

IMPLICIT REAL*8 (A-H,0-Z) 

CCMPLEX*16 GAMN 

REALMS I VEC, ME AN, MHAT, 13, I C , KSDN , KL, KU 
INTEGERS SIZE, SS, HE 

DIMENSION X ( 820 ) , S ( 409 ) , I WK ( 249 5 ) , I VEC (409) 

DATA PI/3. 14159265400/ 

GAMMA 1 = 0 . 0 DO 
G AMMA2 = 0. ODO 
KSDN= 0. ODO 
MEAN= 0. ODO 
V AR I AN= 0 . 0 DO 

NN = SIZE - IDINT (DFL0AT(SIZE)/2) - 1 
SS= NN -1 
HE = SIZE -1 
JE = SIZE - 2 

CALCULATE MEAN AND VARIANCE OF INTERVALS 
DO 10 J = 1 , S I Z E 
MEAN = MEAN + X ( J ) 
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VAR I AN = VARIAN + X(J)**2 
10 CONTINUE 

MEAN = MEAN/SIZE 

VARIAN = ( VARIAN - SIZE * MEAN**2) / (SIZE - 1) 



( X( J+l ) -MEAN ) 
(X( J+2J-MEAN) 



(X( JJ-MEAN) 
( X( J ) -MEAN) 



CALCULATE ESTIMATES OF M AND BETA 
DO 40 J =1 » HE 
GAMMA1 = GAMMA 1 + 

40 CCNTINUE 

DO 50 J = 1 * JE 
GAMMA2 = GAMMA2 + 

50 CONTINUE 

GAMMA1 = GAMMA 1/SIZE 
G AM MA2=GAMMA 2/SIZE 
BHAT = GAMMA2 / GAMMA1 
MHAT = GAMMA 1**2 / GAMMA2 

COMPUTE FINITE FOURIER TRANSFORM OF INTERVAL DATA 
CALL FFTR(X, GAMN,SIZE, IWK) 

CALCULATE PER I ODOG RAM 
DO 20 J =3 » SIZE, 2 
I = ( J-l ) /2 

IVEC( I )=(X( J)**2 + X( J + l )**2) /(PI * SIZE) 

20 CONTINUE 

TEST PERI ODOGRAM FIT TO ESTIMATED POWER SPECTRAL DENSITY 
OMEGA = DCGS ( 2 * PI /SIZE) 

PHAT= (VARI AN+2*MHAT*BHAT* (OMEGA- BHAT) /( l+BHAT**2-2* 
CBF.AT *0M EG A ) ) / P I 
S(l) = IVEC(1)/PHAT 
DO 60 J =2 , NN 

OMEGA = DCOS ( 2 * PI * J / SIZE) 

PHAT=(VARIAN+2*MHAT*BHAT*(0MEGA-BHAT )/ ( 1 +BFAT**2-2* 
CBFAT*0MEGA ) ) /PI 
S(J)= S(J-l) + I V EC ( J ) / PHAT 
60 CONTINUE 

DC 70 J=1,SS 
S ( J ) = S(J)/S(NN) 

KL= DABS(S(J)-(DFLOAT( J-1J/NN) ) 

KU= DAoS( S ( J )-( D FLOAT ( J )/NN ) ) 

KSDN = DMAX1(KSDN,KL,KU) 

70 CONTINUE 

KSDN=KS DN*DS QRT (DFLOAT (NN) ) 

80 RETURN 
END 



SUBROUTINE ESTIM8 

PURPOSE: 

THIS SUBROUTINE USES THE METHOD OF SCORING TO 
STABILIZE ESTIMATES OF THE PARAMETERS FOR A 
UNIVARIATE TwO-STATE SEMI-MARKOV MODEL. 

USAGE: 

1. ARGUMENTS: 



M1,M2,S2,A1,A2 - 



INPUT INITIAL ESTIMATES FOR 
THE MEAN OF TYPE 1 INTERVALS, 
MEAN AND STANDARD CEVIATICN 
OF TYPE 2 INTERVALS AND THE 
TRANSITION PROBABILITIES - 
ALL REAL*8 



S IZE - INPUT NUMBER 
(INTEGER) 



OF INTERVAL DATA POINTS 



33 



ooooonoooooooononooo 



on no on nooooooooooooooooooo 



MEAN - INPUT MEAN OF INTERVAL DATA (REAL*8) 

I VEC - INPUT VECTOR OF P ER I ODOGRAM POINTS 
(REALMS) 

ITER8 - INPUT NUMBER OF ITERATIONS DESIRED 

PRIOR TO TERMINATION IF NO CONVERGENCE 
( INTEGER) 

CONVRG - CONVERGENCE CRITERION FOR TERMINATION 
( RE AL*8 ) 10 . 0E- 10 RECOMMENDED 

2. SUBROUTINES REQUIRED: 

DMINV (IBM ROUTINE) 

DTERM (IBM ROUTINE) 

3. IF NO. SCALING IS DESIRED, SET MEAN=1.0DO . 



SUBROUT INE EST IM8 ( M 1 , M2 , S2 , A1 , A 2 , S IZ E , M E AN , I VEC , I TER8 , 
CCCNVRG ) 

IMPLICIT REAL*8 (A-Z) 

INTEGERS SIZE,NN,J,NE,JE,HE, I ,L,M,N, ITE , ITER8 
DIMENSION AVEC ( 5 ) ,AMAT(5,5) , MV EC ( 5 ) , MMAT (5,5),BVEC(5) 
DIMENSION BMAT(5*5),L(5)»LVEC(5),LMAT(5,5) » CELT ( 5 ) 
DIMENSION LLMAT (5,5) , M( 5) , PVEC ( 5 ) , IV EC (409) 

DATA PI /3. 141592654D0/, N/5/ 

I T E=0 

NN=SIZE-IDINT (DFLOAT( SIZE)/2)-l 

ZERC-OUT VECTORS AND MATRICES 
90 DO 100 J=1 , 5 

AVEC(J) = O.ODO 
BVEC(J) = O.ODO 
LVEC(J) = O.ODO 
MVEC(J) = O.ODO 
DO 110 1=1,5 
AMAT ( I , J) = O.ODO 
BMAT ( I, J) = O.ODO 
LMAT ( I , J ) = O.ODO 

L LMAT ( I , J ) = O.ODO 
MMAT ( I , J) = O.ODO 
110 CONTINUE 
100 CONTINUE 

COUNT ITERATIONS 
ITE = ITE + 1 
WRITE (6,345) ITE 

CALCULATE ELEMENTS OF EQUATIONS 
CC1 = 1 - A1 
CC2 = 1 - A2 
CC3 = CC1 + CC2 
CC4 = ( Ml**2) - Ml 
CC5 = Ml - M2 
CC6 = CC1 * CC2 

AE = ( ( CC2 * CC4 ) + (CC1 * S2**2))/CC3 
AV EC ( 1 ) = MEAN *< ( 1-2*M1 ) *(-CC2 ) ) /CC3 
AVEC ( 3 ) = MEAN * (2* S2 * CC1J/CC3 
A VEC ( 4 ) = ( CC2 * ( CC 4 - S2**2))/ (CC3**2) 

A VEC ( 5 ) = ( (-CC1)#(CC4 -S2**2))/ (CC3**2) 

AMAT (1,1)= MEAN *(2 * CC2J/CC3 
AMAT (1,4)= AVEC(l) / ( CC3 * MEAN) 

AMAT ( 1 , 5 ) = (CC1 *(1-2*M1) )/(CC3**2) 

AMAT (3,3)= MEAN =M2* CC1)/ CC3 

AM AT ( 3 , 4 ) = (2 * S2 * (-CC2 ) ) / CC3 #*2 

A M A T ( 3 , 5 ) = A V EC ( 3 ) / ( CC3 * MEAN) 

AM AT ( 4 , 4 ) = (2 * AVEC(4)) / CC3 
AM AT ( 4 , 5 ) = ((A1-A2) * ( CC4 - S2**2))/ CC3#*3 
AM AT( 5 , 5 ) = (2 * AVEC(5))/CC3 
ME = (CC6 * CC5**2 )/CC3**2 
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MVEC( 1 ) = MEAN * (2 * CC5 * CC6 ) / CC3**2 
M V EC ( 2 ) = -MVEC( 1 ) 

MVEC(4) =(CC2*(A2-A1)*(CC5**2) )/CC3**3 
MVEC(5) =(CC1*(A1-A2)*(CC5**2) ) /CC3**3 
MMAT(1 ,1) = MEAN*(2*CC6)/CC3**2 
MMAT( 1, 2)=-MMAT( 1,1) 

MMAT( 1,4)= (2*CC2*(A2-A1)*CC5)/CC3**3 
MMAT (2,4) =-MMAT ( 1,4) 

MM AT ( 1, 5 ) = (2*CC1*(A1-A2)*CC5)/CC3**3 
MM AT (2,5)=- MM AT ( 1,5) 

MMAT(2,2)=MEAN*(2*CC6)/CC3**2 

MMAT (4, 4)= ( ( (6*A2)+< < -4 )* A2**2 ) + ( 2*A1*A2 )+( (-2)*Al)-2) 
C* (CC5**2) ) /CC3**4 

MMAT ( 4, 5 ) = ( ( (2*( 1-A1-A2+2*A 1*A2 ) )-Al**2-A2**2 ) *CC5**2 ) 
C/CC3**4 

MMAT( 5,5) =( ( (6*Al)-4*Al**2+2*Al*A2-2*A2-2 ) *(005**2 ) ) 
0/003**4 
DO 120 J=1 ,3 
DO 130 1=1,5 

AMAT ( J , I ) = AMAT ( J , I ) *MEAN 
MMAT ( J , I ) = MMAT ( J , I )*MEAN 
130 CCNTINUE 
120 CCNTINUE 

BETA=A1+A2-1.0DO 
L IKE= 0 . ODO 
DC 140 N E = 1 , N N 

COSINE = DOGS ( 2 * PI * NE/ SIZE) 

CC7 = (1 + BET A**2 ) - 2 * BETA * COSINE 
BE = 1 + ( 2 *( (BETA * COSI NE) -BETA**2 ) J/CC7 
B V EC ( 4 ) = 2*( ( ( 1+BETA**2)*CC SINE) -2*BETA) /CC 7**2 
B VEC( 5 ) =BV EC ( 4 ) 

BMAT (4,4) = 2* ( 6*B ET A** 2-2* CCS I N E*BET A**3 + 4*CGS I NE**2-6* 
CBETA*C0SINE-2)/CC7**3 
BMAT (4,5) = BMAT (4,4) 

EMAT ( 5 , 5 ) = BMAT (4,4) 

CALCULATE ESTIMATE OF POWER SPECTRAL DENSITY 
PE= (AE + ME *BE ) / PI 
CC8 = ( I VEC (NE ) -PE )/ PE**2 
CC9 = (PE - 2 * IVEC(NEJ) / PE**3 

CALCULATE FIRST-ORDER PART I ALS OF PSD 
DO 150 J = 1 , 5 

PVEC(J) = (AVEC(J) + MVEC(J )*BE + BVEC(J)* ME) / PI 
150 CONTINUE 

CALCULATE LIKELIHOOD FUNCTION VALUE 

LIKE = LIKE - (IVEC(NE) / PE ) - DLOG(PE) 

CALCULATE SCORES 
DO 160 JE = 1 , 5 

LVEC(JE)=LVEC (JE) + CC8 * PVEC(JE) 

DO 170 H E = J E , 5 

CALCULATE SECOND-ORDER PARTI ALS OF PSD 

PJH=AMAT(JE,HE) + MMAT ( JE, HE ) * BE + MVEC(JE) * 
CBVEC(HE)+MVEC ( HE ) *B VEC ( J E ) + BMAT ( JE , HE ) * ME 

CALCULATE SECOND-ORDER PART I ALS OF LIKELIHGGD FUNCTION 
LMAT ( J E , HE ) = LMAT (JE,HE )-CC 9*PV EC ( J E ) *P VEC (HEJ-CC8* 
CPJh/PI 

CALCULATE EXPECTED VALUE OF LMAT 

LLMAT ( J E , HE) =LLMAT ( JE , HE ) +PVEC ( J E ) *PV EC ( HE ) /PE**2 
170 CONTINUE 
160 CCNTINUE 
140 CGNTINUE 

FILL IN LOWER TRIANGLE OF MATRICES 
DO 180 JE = 1 , 4 
J = J E + 1 
DG 190 HE = J , 5 
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LMAKHE , JE )=LMAT{ JE,HE ) 
LLMAKHE, J E ) =LLM AT { JE » HE ) 

190 CONTINUE 
180 CONTINUE 

WRI TE ( 6 , 360 ) 

DO 200 1=1,5 

WRITE(6,370) ( LM AT ( I , J ) , J= 1 , 5 ) 

DO 185 J=1 , 5 

AMAT ( I , J J =-LMAT ( I, J) 

185 CONTINUE 
200 CONTINUE 

CALCULATE INVERSE AND DETERMINANT 
CALL DM INV ( LMAT ,N,D,L,M) 
WRITE(6,380) D 



OF LMAT 



210 



WRITE ( 6 , 390 ) 
DO 210 1=1.5 
W R ITE ( 6 , 370 ) 
CONTINUE 



( LM AT ( I , J ) ,J=1, 5) 



CALCULATE ADDITIVE 
DO 240 1=1,5 
D ELT ( I ) = 0.0D0 
DO 250 J=1 , 5 
DELTII ) = DELT ( I ) 
250 CONTINUE 
240 CONTINUE 



INCREMENTS TO ESTIMATES 



+ LM AT ( I , J ) * LVEC(J) 



INCREMENT ESTIMATES 



M 1= Ml + 

M2= M2 + 

S2= S2 + 

A 1= A 1 + 

A2= A2 + 
WRITE(6,430) 
W R I TE ( 6 , 44 0 ) 
WRITE ( 6 ,450) 



DELT ( 1 ) 
D ELT ( 2 ) 
DELT (3) 
DELT ( 4) 
DELT ( 5 ) 



* 

* 



MEAN 

MEAN 

MEAN 



LIKE 

( L VEC ( J ) , J = 1 , 5) 
Ml , M2 , S2 , A 1 , A2 



TEST FOR CONVERGENCE 

I F ( DMAX 1 ( DABS (DELT ( 1 ) ) , DABS ( DELT ( 2 ) ) , DA BS ( CELT ( 3 ) ) , 
CDABS(DELT (4) ) ,DABS(DELT(5) ) ) .LE.CCNVRG) GO TO 260 
I F ( DMAX 1 ( DABS (LV EC ( 1 ) ) , DABS ( LV EC ( 2 ) ) , DAB S ( L VEC ( 3 ) ) , 
CDABS(LVEC(4) ) ,DABS(LVEC(5) ) ) . L E . CONVRG ) GO TO 260 

TEST FOR NUMBER OF ITERATIONS 
255 I F ( ITE.LT. ITER8) GO TO 270 
280 WRITE (6,455) CONVRG 
GO TO 265 
270 WRI TE ( 6 ,460) 

CONTINUE IF NECESSARY 
GO TO 90 
260 WRITE(6 ,470) 

TEST LMAT FOR NEGATIVE DEFINITENESS 
265 DC 275 1=1,5 
DO 267 JE=1 , 5 
DC 266 HE= 1 , 5 
BMAT(JE,HE)=AMAT ( JE, HE ) 

CONTINUE 
CONTINUE 

CALL DTERMU , BMAT , D, N ) 

A V EC ( I ) =D 
CONTINUE 

I F ( AVEC ( 1 ) . LT.O.ODO. AND. AVEC(2 ) . GT . 0 . ODO . ANC . AV EC ( 3 ) 
C.LT .O.GDO.AND.AVEC{4) .GT. 0. ODO . AND.A VEC ( 5 ) . LT . O.ODO) 
C GC TO 290 
285 WRITE(6,480) 

DC 295 1=1,5 
WRITE(6,490) I , AVEC ( I ) 

CCNT INUE 
GO TO 600 



266 

267 



275 



295 
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no on 



290 WRITE(6,500) 
600 WRITE(6,400) 
DO 220 1 = 1,5 
WRITE(6,370) 
220 CONTINUE 



(LLMAT< I, J),J=1,5) 



INVERT EXPECTED VALUE MATRIX 
CALL DMINV(LLMAT ,N,D,L,M) 

WRITE(6,410) D 
WRITEI6 ,420) 

TEST INVERSE FOR POSITIVE DEFINITENESS 
DC 230 1=1,5 

WRITE ! 6 , 370 ) ( LLMAT ( I , J ) , J = 1 , 5) 

DO 225 J = 1 » 5 

AMAT( I, J)=LLMAT(I , J) 

225 CONTINUE 
230 CONTINUE 

DO 610 1=1,5 

DO 620 JE= 1, 5 

DC 630 HE=1 , 5 

BMAT( JE , HE )=AMAT ( JE , HE ) 

630 CONTINUE 
620 CONTINUE 

CALL DTERMU ,BMAT ,D,N) 

AVEC( I ) =D 
610 CONTINUE 

I F( DMIN 1 ( AVEC( 1 ) , AVEC(2) , AVEC ( 3 ) , A VEC ( 4 ) , A VEC ( 5) ) .GT . 
CO.ODO) GO TO 660 



640 



650 

660 

345 

360 

370 

380 

390 

400 

410 



WRITE ( 6, 510) 

DG 650 1=1,5 

WRI TE ( 6 , 520) I,AVEC(I) 

CONT INUE 
RETURN 

WRITE(6,530) 

RETURN 

FORMAT ( 1H1 ,//////////, • ITERATION NUMBER', 14) 
FORMAT!//, • NEGATIVE MATRIX OF SECOND P ART I ALS ' 
FORMAT (/,5D20.10) 



/) 



DETERMINANT OF MATRIX' 
INVERSE MATRIX' ,/) 
INFORMATION MATRIX',/) 
DETERMINANT OF INFORMATION 



' INVERSE INFORMATION MATRIX',/) 



FORMAT! //, 

FORMAT!//, 

FORMAT ( // , 

FORMAT! //, 

CD20.10) 

420 FORMAT!//, 

430 FORMAT!//, 

440 FORMAT!//, 

450 FORMAT! //, 

455 FORMAT!//, 

C' REACHED 
CD18.10,' ) 

460 FORMAT!//, 

470 FORMAT!//, 

480 FORMAT! 1H1, /////, • MATRIX 
C ' NEGATI VE DEF INI TE' ,/ ) 

490 FORMAT! /,' DETERMINANT CF 
C'CF ORDER', 12,' OF MATRIX 
CD20.10) 

500 FORMAT! 1H1, /////, • MATRIX 
C' NEGATI VE DEFINITE' ) 



,//, D20 .10) 



MATRIX' ,// 



LOG LIKELIHOOD FUNCTION 
' SCORES =' ,5020. 10) 

' PARAMETERS =',5D20.10) 

' NUMBER OF ITERATIONS REQUESTED 
AMD',//,' CONVERGENCE CRITERION 
WAS NOT ACHIEVED. ' , // 

' CONTINUING') 



VALUE = ' , D20 .10) 



RUN 



WAS' , 

i ' , 

TERMINATED' ) 



COVERGENCE' 

i 



) 



OF SECOND PART I ALS IS NOT ' 



PRINCI PAL 
OF SECOND 



MINOR SUBMATRIX 
PART IALS = ' , 



OF SECOND PART IALS IS ' 



510 FORMAT! ///,' INVERSE INFORMATION MATRIX IS NOT ' 
C' POSITIVE DEFINITE' ,/ ) 



52C FORMAT!/, • DETERMINANT OF PRINCIPAL MINOR SUBMATRIX • 
C ' OF ORDER', 12,' OF INVERSE INFORMATION MATRIX = ', 
CD20.10) 

530 FORMAT!///,' INVERSE INFORMATION MATRIX IS POSITIVE ' 
C' CEFINITE' ) 

END 
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